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Abstract. The gravitational clustering of collisionless particles in an expanding 
universe is modelled using some simple physical ideas. I show that it is possible to 
understand the nonlinear clustering in terms of three well defined regimes: (1) linear 
regime; (2) quasilinear regime which is dominated by scale- invariant radial infall and 
(3) nonlinear regime dominated by nonradial motions and mergers. Modelling each 
of these regimes separately I show how the nonlinear two point correlation function 
can be related to the linear correlation function in hierarchical models. This analysis 
leads to results which are in good agreement with numerical simulations thereby 
providing an explanation for numerical results. Using this model and some simple 
extensions, it is possible to understand the transfer of power from large to small 
scales and the behaviour of higher order correlation functions. The ideas presented 
here will also serve as a powerful analytical tool to investigate nonlinear clustering 
in different models. 



1. Introduction 

Consider a collection N point particles, interacting with each other by 
the Newtonian gravity, in an expanding background characterized by 
a scale factor a(t). What can we say about the time evolution of such 
a system? 

This problem is of considerable interest for several reasons. To begin 
with, the behaviour of large number of particles interacting via Newto- 
nian gravity poses a formidable challenge to the usual methods of statis- 
tical mechanics (T. Padmanabhan, 1990). So, purely from an academic 
point of view, this seems to be a challenging but solvable problem. 

Secondly, this problem might even have some practical interest. 
There is considerable evidence that the universe is dominated by col- 
lisionless non-baryonic dark matter particles. In that case, they will 
play a key role in the formation of large scale structures. If the length 
scales of interest are (i) small compared to Hubble radius but (ii) large 
compared to the scales at which non-gravitational processes are signif- 
icant, then the system of dark matter particles constitutes an example 
in which the question raised in the first paragraph becomes relevant. 
In fact, most of the work in this subject has been inspired by consid- 
erations of structure formation. 
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A brute force method for solving this problem relies on numerical 
simulations. In such an approach one starts with large number of parti- 
cles distributed nearly uniformly and calculates the future trajectories 
by a suitable numerical algorithm. This approach, however, does not 
lead to genuine understanding unless we supplement it with some ana- 
lytic modelling. In this talk, I shall outline how one can make analytic 
progress in the problem of nonlinear gravitational clustering and there- 
by reproduce the key features of numerical simulations (Padmanabhan, 
1996 a,b). 

2. Modelling the gravitational clustering 

Consider a system of particles distributed homogeneously, on the aver- 
age, with a mean density p(t). This uniform density will cause an 
expansion of the universe and the proper distance r = a(i)x between 
particles will increase with time. If the distribution was not strictly 
uniform, then the perturbations in the density will act as local centres 
of clustering. A region with overdensity will accrete matter around it 
while an underdense region will repel matter in its surroundings. As a 
result, perturbations in density will tend to grow and, when the densi- 
ty contrast is of order unity, these cluster centres will exert significant 
influence on the evolution. Particles in a highly overdense regions will 
evolve essentially under their own self-gravity and will tend to form 
gravitationally bound systems. 

When these density perturbations are small, it is possible to study 
their evolution using linear theory. But once the density contrast becomes 
comparable to unity, linear perturbation theory breaks down and one 
must use N-body simulations to study the growth of perturbations. 
While these simulations are of some value in making concrete predic- 
tions for specific models, they do not provide clear physical insight into 
the process of non-linear gravitational dynamics. To obtain such an 
insight into this complex problem, it is necessary to model the gravi- 
tational clustering of collisionless particles using simple physical con- 
cepts. I shall develop one such model in this section, which - in spite of 
extreme simplicity - reproduces the simulation results for hierarchical 
models fairly accurately. Further, this model also provides insight into 
the clustering process and can be modified to take into account more 
complicated situations. 

The paradigm for understanding the clustering is based on the well 
known behaviour of a spherically symmetric overdense region in the 
universe. In the behaviour of such a region, one can identify three dif- 
ferent regimes of interest: (1) In the early stages of the evolution, when 
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the density contrast is small, the evolution is described by linear the- 
ory. (2) Each of the spherical shells with an initial radius X{ can be 
parametersed by a mass contained inside the shell, M(xj), and the 
energy, E{xi) for the particular shell. Each shell will expand to a maxi- 
mum radius x max oc M/\E\ and then turn around and collapse. Such a 
spherical collapse and resulting evolution allows a self similar descrip- 
tion (Filmore & Goldreich, 1984; Bertshinger, 1985) in which each shell 
acts as though it has an effective radius proportional to x max . This will 
be the quasilinear phase. (3) The spherical evolution will break down 
during the later stages due to several reasons. First of all, non radi- 
al motions will arise due to amplification of deviations from spherical 
symmetry. Secondly, the existence of substructure will influence the 
evolution in a non-spherically symmetric way. Finally, in the real uni- 
verse, there will be merging of such clusters [each of which could have 
been centres of spherical overdense regions in the begining] which will 
again destroy the spherical symmetry. This will be the nonlinear phase. 

The description given above is sufficiently well known that one may 
suspect it can not lead to any insight into the problem. In particular, 
structures observed in the real universe are hardly spherical. I will show 
that it is, however, possible to model the above process in a manner 
which allows direct generalisation to the real universe. 

To do this we will begin by studying the evolution of system start- 
ing from a gaussian initial fluctuations with an initial power spectrum, 
Pin(k). The fourier transform of the power spectrum defines the corre- 
lation function £(a, x) where a oc t 2 / 3 is the expansion factor in a uni- 
verse with fi = 1. It is more convenient for our purpose to work with 
the average correlation function inside a sphere of radius x, defined by 

3 f x 

€(a,x) = -3 / £(a,y)y 2 dy (1) 
x^ Jo 

In the linear regime we have £l(ci, x) oc a 2 £j n (aj, x). In the quasilinear 
and nonlinear regimes, we would like to have prescription which relates 
the exact ^ to the mean correlation function calculated from the linear 
theory. One might have naively imagined that £(a, x) should be related 
to x). But one can convince oneself that the relationship is likely 
to be nonlocal by the following analysis: 

Recall that, the conservation of pairs of particles, gives an exact 
equation satisfied by the correlation function (Peebles, 1980): 

where v(t, x) denotes the mean relative velocity of pairs at separation x 
and epoch t. Using the mean correlation function £ and a dimensionless 
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pair velocity h{a,x) = — (y/ax), equation (g) can be written as 
o In a a In x 

This equation can be simplified by first introducing the variables 

A = In a, X = lnx, D(X,A) = ln(l +£) (4) 
in terms of which we have (Nityananda and Padmanabhan, 1994) 
dD dD 

— - h (A,X)—=3h(A,X) (5) 

Introducing further a variable F = D + 3X, (|5|) can be written in a 
remarkably simple form as 

dF dF 

— - h(A,X)— = (6) 

The charecteristic curves to this equation - on which F is a constant - 
are determined by (dX/dA) = —h(X, A) which can be integrated if h 
is known. But note that the charecteristics satisfy the condition 

F = 3X + D = ln[z 3 (l + £)] = constant (7) 

or, equivalently, 

x 3 (l + = l 3 (8) 

where I is another length scale. When the evolution is linear at all the 
relevant scales, ^ < 1 and I « x. As clustering develops, £ increases and 
x becomes considerable smaller than /. It is clear that the behaviour 
of clustering at some scale x is determined by the original linear power 
spectrum at the scale / through the "flow of information" along the 
charesteristics. This suggests that we should actually try to express 
the true correlation function £(a,x) in terms of the linear correlation 
function evaluated at a different point. 

Let us see how we can do this starting from the quasilinear regime. 
Consider a region surrounding a density peak in the linear stage, around 
which we expect the clustering to take place. It is well known that 
density profile around this peak can be described by 

p(x)tt Pbg [l + Z(x)} (9) 

Hence the initial mean density contrast scales with the initial shell 
radius / as 5i(l) oc £l(0 in the initial epoch, when linear theory is valid. 
This shell will expand to a maximum radius of x max oc l/5i oc Z/£l(0- ^ n 
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scale-invariant, radial collapse, models each shell may be approximated 
as contributing with a effective radius which is propotional to x max . 
Taking the final effective radius x as proportional to x max , the final 
mean correlation function will be 

M I 3 

Zq L (x) oc p oc — oc (p/ ^ (0)3 oc l L (lf (10) 

That is, the final correlation function £ql at x is the cube of initial 
correlation function at I where I 3 oc x 3 ^\ oc x 3 ^ql(x). This is in the 
form demanded by (|8|) if £ 3> 1. Note that we did not assume that the 
initial power spectrum is a power law to get this result. 

In case the initial power spectrum is a power law, with c^l oc x-(™+ 3 ) ; 
then we immediately find that 

£ Q L OC x -3(n+3)/(n+4) (n) 

[If the correlation function in linear theory has the powerlaw form 
£l oc x~ a then the process desribed above changes the index from 
a to 3a/(l + a). We shall comment more about this aspect later]. 
For the power law case, the same result can be obtained by more 
explicit means. For example, in power law models the energy of spher- 
ical shell with mean density 5(xi) oc x^ will scale with its radius as 
E oc G5M(xi)/xi oc G5xf oc x\~ h . Since M oc xf, it follows that the 
maximum radius reached by the shell scales as x max oc (M/E) oc x\ +b . 
Taking the effective radius as x = x e ff oc x\ +h , the final density scales 
as 

poc J oc -J^ oc sr* oc s-3W+» (12) 



In this quasilinear regime, £ will scale like the density and we get 
|ql oc a;-3b/(i+6)_ fp^g i nc ^ ex b can kg related to n by assuming the the 
evolution starts at a moment when linear theory is valid. Since the grav- 
itational potential energy [or the kinetic energy] scales as E oc x i ( n+1 ^ 
in the linear theory, it follows that b = n + 3. This leads to the corre- 
lation function in the quasilinear regime, given by (|ll|) . 

If O = 1 and the initial spectrum is a power law, then there is no 
intrinsic scale in the problem. It follows that the evolution has to be 
self similar and £ can only depend on the combination q = xa~ 2 ^ n+3 \ 
This allows to determine the a dependence of ^ql by substituting q for 
x in (jll|). We find 

£ QL {a, x) oc a 6 /(" +4 )x- 3 ( n+3 )/( n+4 ) (13) 
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We know that, in the linear regime, £ = £l oc a 2 . Equation (|ll|) shows 
that, in the quasilinear regime, £ = £ql oc a 6 ^ n+4 \ Spectra with n < 
— 1 grow faster than a 2 , spectra with n > — 1 grow slower than a 2 and 
n = — 1 spectrum grows as a 2 . 
Direct algebra shows that 

£ QL (a, x) oc [h(a,l)f (14) 

reconfirming the local dependence in a and nonlocal dependence in 
spatial coordinate. This result has no trace of original assumptions 
[spherical evolution, scale-invariant spectrum ....] left in it and hence 
once would strongly suspect that it will have far general validity. 

Let us now proceed to the third and nonlinear regime. If we ignore 
the effect of mergers, then it seems reasonable that virialised systems 
should maintain their densities and sizes in proper coordinates, i.e. 
the clustering should be "stable". This would require the correlation 
function to have the form £jvx(o,x) = a 3 F(ax). [The factor a 3 arising 
from the decrease in background density]. Prom our previous analysis 
we expect this to be a function of £j,(a, Z) where I 3 « x 3 ^nl{cl,x). Let 
us write this relation as 

tNLia, x) = a 3 F(ax) = U[£ L (a, I)] (15) 

where U[z] is an unknown function of its argument which needs to be 
determined. Since linear correlation function evolves as a 2 we know 
that we can write £l(o> I) = a 2 Q[l 3 ] where Q is some known function 
of its argument. [We are using I 3 rather than I in defining this function 
just for future convenience of notation]. In our case I 3 = x 3 £nl(cl, x) = 
(ax) 3 F(ax) = r 3 F(r) where we have changed variables from (o, x) to 



(a, r) with r = ax. Equation (15) now reads 



a 3 F(r) = U[£ L (a, I)] = U[a 2 Q[l 3 ]] = U[a 2 Q[r 3 F(r)]] (16) 

Consider this relation as a function of a at constant r. Clearly we need 
to satisfy £/[cia 2 ] = C2a 3 where c\ and ci are constants. Hence we must 
have 

U[z] oc z 3/2 . (17) 
Thus in the extreme nonlinear end we should have 

6vL(a,x)oc[a(a,Z)] 3 / 2 (18) 

[Another way deriving this result is to note that if £ = a 3 F(ax), then 
h = 1. Integrating (||) with appropriate boundary condition leads to 
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(18) .] Once again we did not need to invoke the assumption that the 



spectrum is a power law. If it is a power law, then we get, 

U(a, x) oc a^x^- 7 = *j2±§. (19 ) 

(n + 5) 

This result is based on the assumption of "stable clustering" and 
was originally derived by Peebles (Peebles, 1965). It can be directly 
verified that the right hand side of this equation can be expressed in 
terms of q alone, as we would have expected. 

Putting all our results together, we find that the nonlinear mean 
correlation function can be expressed in terms of the linear mean cor- 
relation function by the relation: 

(h{a,l) (fora<l, e<i) 

f (o, x) = < 3 ( for 1 < & < 5 - 85 > 1 < £ < 20 °) ( 20 ) 

{ 14.14f L (a, if' 2 (for 5.85 < f £ , 200 < f) 

The numerical coefficients have been determined by continuity argu- 
ments. We have assumed the linear result to be valid upto £ = 1 and 
the virialisation to occur at £ ~ 200 which is result arising from the 
spherical model. The exact values of the numerical coefficients can be 
obtained only from simulations. 

The true test of such a model, of course, is N-body simulations and 
remarkably enough, simulations are very well represented by relations 
of the above form. The simulation data for CDM, for example, is well 
fitted by (Padmanabhan etal., 1996): 

fa(a,0 (for £ L < 1.2, £< 1.2) 

£(a, x) = I £ L (a, if (for 1< £ L < 5, 1< £ < 125) (21) 

{ 11.7^(0, if 2 (for 5 < 125 < £) 

which is fairly close to the theoretical prediction. [The fact that numer- 
ical simulations show a correlation between £(a, x) and £i(a, I) was 
originally pointed out by Hamilton et al. (1991) who, however, tried to 
give a multiparameter fit to the data. This fit has somewhat obscured 
the simple physical interpretation of the result though has the virtue 
of being very accurate for numerical work.] 

A comparison of ( |20| ) and ( |2l|) shows that the physical processes 
which operate at different scales are well represented by our model. 
In other words, the processes descibed in the quasilinear and nonlinear 
regimes for an individual lump still models the average behaviour of the 
universe in a statistical sense. It must be emphasised that the key point 
is the "flow of information" from / to x which is an exact result. Only 
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when the results of the specific model are recast in terms of suitably 
chosen variables, we get a relation which is of general validity. It would 
have been, for example, incorrect to use spherical model to obtain rela- 
tion between linear and nonlinear densities at the same location or to 
model the function h. 

It may be noted that to obtain the result in the nonlinear regime, 
we needed to invoke the assumption of stable clustering which has not 
been deduced from any fundamental considerations. In case mergers 
of structures are important, one would consider this assumption to be 
suspect (see Padmanabhan et al., 1996). We can, however, generalise 
the above argument in the following manner: If the virialised systems 
have reached stationarity in the statistical sense, the function h - which 
is the ratio between two velocities - should reach some constant value. In 
that case, one can integrate (g) and obatin the result £nl = a 3h F(a h x). 
A similar argument will now show that 

t NL (a,x)<x[t L (a,l)? h/2 (22) 

in the general case. For the power law spectra, one would get 

Simulations are not accurate enough to fix the value of h; in particular, 
the asymptotic value of h could depend on n within the accuracy of 
the simulations. It may be possible to determine this dependence by 
modelling mergers in some simplified form. 

If h = 1 asymptotically, the correlation function in the extreme 
nonlinear end depends on the linear index n. One may feel that physics 
at highly nonlinear end should be independent of the linear spectral 
index n. This will be the case if the asymptotic value of h satisfies the 
scaling 

h = (24) 
n + 3 

in the nonlinear end with some constant c. Only high resolution numer- 
ical simulations can test this conjecture that h{n + 3) = constant. 

It is possible to obtain similar relations between £(a, x) and £i,(a, I) 
in two dimensions as well. In 2-D the scaling relations turn out to be 

{^l(o, I) (Linear) 
£l(°>0 2 (Quasi-linear) (25) 
£l(o.,1) (Nonlinear) 

For power law spectrum the nonlinear correction function will ^jvl(O) x ) - 
a^x-i with 7 = 2(n + 2) /(n + 4). 
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If we generalize the concept of stable clustering to mean constancy 
of h in the nonlinear epoch, then the correlation function will behave 
as £atl(o, x) = a 2h F(a h x). In this case, if the spectrum is a power law 
then the nonlinear and linear indices are related to 

_ 2h{n + 2) 
7 - 2 + h(n + 2) [ ] 

All the features discussed in the case of 3 dimensions are present here 
as well. For example, if the asymptotic value of h scales with n such 
that h(n + 2) = constant then the nonlinear index will be independent 
of the linear index. (Numerically it would be lot easier to test this result 
in 2-D rather than in 3-D; work is in progress to test these results). 

We shall now consider some applications and further generalisations 
of our model. 



3. Critical Index and power transfer 

Given a model for the evolution of the power spectra in the quasilinear 
and nonlinear regimes, one could explore whether evolution of gravi- 
tational clustering possesses any universal charecteristics. For example 
one could ask whether a complicated initial power spectrum will be 
driven to any particular form of power spectrum in the late stages of 
the evolution. 

One suspects that such a possibility might arise because of the fol- 
lowing reason: We saw in the last section that [in the quasilinear regime] 
spectra with n < — 1 grow faster than a 2 while spectra with n > — 1 
grow slower than a 2 . This feature could drive the spectral index to 
n = n c ~ — 1 in the quasilinear regime irrespective of the initial index. 
Similarly, the index in the nonlinear regime could be driven to n ~ —2 
during the late time evolution. So the spectral indices —1 and —2 are 
some kind of "fixed points" in the quasilinear and nonlinear regimes. 
Speculating along these lines, we would expect the gravitational cluster- 
ing to lead to a "universal" profile which scales as x^ 1 at the nonlinear 
end changing over to x~ 2 in the quasilinear regime. 

This effect can be understood better by studying the "effective" 
index for the power spectra at different stages of the evolution. These 
are plotted in figure 1. The three panels of figure 1 illustrate features 
related to the existence of fixed points in a clear manner. In the top 
panel we have plotted index of growth n a = (91n£(a, x)/dlna) x as 
a function of £ in the quasilinear regime obtained from our scaling 
relations. Curves correspond to an input spectrum with index n = 
—2, —1, 1. The dashed horizontal line at n a = 2 represents the linear 
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growth rate. An index above this dashed horizontal line will represent 
a rate of growth faster than linear growth rate and the one below will 
represent a rate which is slower than the linear rate. It is clear that 
- in the quasilinear regime - the curve for n = — 1 closely follows the 
linear growth while n = — 2 grows faster and n = 1 grows slower; so the 
critical index is n c « — 1. The curves are based on the fitting formula 
due to Hamilton et al, 1991. 

The second panel of figure 1 shows the effective index func- 
tion of the index n of the original linear spectrum at different levels of 
nonlinearity labelled by £ = 1, 5, 10, 50, 100. We see that in the quasi- 
linear regime, n a > 2 for n < — 1 and n a < 2 for n > — 1. 

The lower panel of figure 1 shows the slope n x = — 3— (<91n£/<91nx) a 
of £ for different power law spectra. It is clear that n x crowds around 
n c ~ —1 in the quasilinear regime. If perturbations grow by gravitaional 
instability, starting from an epoch at which initial < 1 at all scales, 
then it can be shown that n x at any epoch must satisfy the inequality 

n x < (3/f) . (27) 

This bounding curve is shown by a dotted line in the figure. This pow- 
erful inequality shows that regions of strong nonlinearity [with f » 1] 
should have effective index which is close to or less than zero. 

The index n c = — 1 corresponds to the isothermal profile with £(a, x) = 
a 2 x~ 2 and has two interesting features to recommend it as a candidate 
for fixed point: 

(i) For n = — 1 spectra each logarithmic scale contributes the same 
amount of correlation potential energy. If the regime is modelled by 
scale invariant radial flows, then the kinetic energy will scale in the same 
way. It is conceivable that flow of power leads to such an equipartition 
state as a fixed point though it is difficult prove such a result in any 
generality. 

(ii) It can be shown that scale invariant spherical collapse will change 
the density profile x~ b with an index b to another profile with index 
36/(1 + b). Such a mapping has a nontrivial fixed point for b = 2 
corresponding to the isothermal profile and an index of power specrum 
n = — 1 (see Padmanabhan, 1996a). 

These considerations also allow us to predict the nature of power 
transfer in gravitational clustering. Suppose that, initially, the power 
spectrum was sharply peaked at some scale ko = 2-k/Lq and has a small 
width Ak. When the peak amplitude of the spectrum is far less than 
unity, the evolution will be described by linear theory and there will be 
no flow of power to other scales. But once the peak approaches a value 
close to unity, power will be generated at other scales due to nonlinear 
couplings even though the amplitude of perturbations in these scales 



icgc.tex; 1/02/2008; 3:16; no v.; p. 10 



Nonlinear Gravitational Clustering in Expanding Universe 



11 




- \ 

- . \ 




1 




• \\ 
■ \ ■ 




_ _ _ _ 5 




: n W 




10 




■. '\\ 




50 








_ . . . 100 




















3 -2 


-1 


12 3 






Figure 1. The top panel shows exponent of rate of growth of density fluctuations as 
a function of amplitude. We have plotted the rate of growth for three scale invariant 
spectra n = —2, — 1, 1. The dashed horizontal line indicates the exponent for linear 
growth. For the range 1 < 8 < 100, the n = — 1 spectrum grows as in linear theory; 
n < — 1 grows faster and n > — 1 grows slower. The second panel shows exponent of 
rate of growth as a function of linear index of the power spectrum for different values 
of £ (1,5,10,50,100). These are represented by thick, dashed, dot-dashed, dotted 
and the dot-dot-dashed lines respectively. It is clear that spectra with n Un < — 1 
grow faster than the rate of growth in linear regime and nun > — 1 grow slower. The 
lower panel shows the evolution of index n x = — 3 — (<91n£/<91n:r) a with £. Indices 
vary from n = —2.5 to n = 4.0 in steps of 0.5. The tendency for n x to crowd around 
n c = — 1 is apparent in the quasilinear regime. The dashed curve is a bounding curve 
for the index (n x < 3/|) if perturbations grow via gravitational instability. 
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are less than unity. Mathematically, this can be understood from the 
evolution equations for the density contrast - written in fourier space - 
as : 

4 + 2-4 = ^Gp5 k + Q (28) 
a 



where 5 k (t) is the fourier transform of the density contrast, p is the 
background density and Q is a nonlocal, nonlinear function which cou- 
ples the mode k to all other modes k' (Peebles, 1980). Coupling between 
different modes is significant in two cases. The obvious case is one with 
5 k > 1. A more interesting possibility arises for modes with no initial 
power [or exponentially small power]. In this case nonlinear coupling 
provides the only driving terms, represented by Q in equation (p8|) . 
These generate power at the scale k through mode-coupling, provided 
power exists at some other scale. Note that the growth of power at the 
scale k will now be governed purely by nonlinear effects even though 
Sk « 1. 

Physically, this arises along the following lines: If the initial spectrum 
is sharply peaked at some scale Lq, first structures to form are voids 
with a typical diameter Lq. Formation and fragmentation of sheets 
bounding the voids lead to generation of power at scales L < Lq. First 
bound structures will then form at the mass scale corresponding to Lq. 
In such a model, £n n at L < Lq is nearly constant with an effective index 
of n rj —3. Assuming we can use equation (^) with the local index 
in this case, we expect the power to grow very rapidly as compared to 
the linear rate of a?. [The rate of growth is a 6 for n = —3 and a 4 for 
n = —2.5.] Different rate of growth for regions with different local index 
will lead to steepening of the power spectrum and an eventual slowing 
down of the rate of growth. In this process, which is the dominant one, 
the power transfer is mostly from large scales to small scales. [There is 
also a generation of the k A tail at large scales which we shall not discuss 
here; see Bagla and Padmanabhan, 1996]. 

From our previous discussion, we would have expected such an evo- 
lution to lead to a "universal" power spectrum with some critical index 
n c ~ — 1 for which the rate of growth is that of linear theory - viz., a 2 . 
In fact, the same results should hold even when there exists small scale 
power; recent numerical simulations dramatically confirm this predic- 
tion and show that - in the quasilinear regime, with 1 < 5 < 100 - power 
spectrum indeed has a universal slope [see Bagla and Padmanabhan, 
1996]. 
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4. Further generalizations 

The ideas presented here can be generalised in two obvious directions 
(see Munshi and Padmanabhan, 1996): (i) By considering peaks of dif- 
ferent heights, drawn from an initial gaussian random field, and aver- 
aging over the probability distribution one can obtain a more precise 
scaling relation, (ii) By using a generalised ansatz for higher order cor- 
relation functions, one can attempt to compute the Sn parameters in 
the quasilinear and nonlinear regimes. I shall briefly comment on the 
results of these two generalisations. 

(i) The basic idea behind the model used in section 2 can be described 
as follows: Consider the evolution of density perturbations starting 
from an initial configuration, which is taken to be a realisation of a 
Gaussian random field with variance a. A region with initial density 
contrast <5j will expand to a maximum radius Xf = Xi/Si and will con- 
tribute to the two-point correlation function an amount proportional 
to (xi/xf) 3 = Sf. The initial density contrast within a randomly placed 
sphere of radius X{ will be vo(xi) with a probability proportional to 
exp(— v 2 /2). On the other hand, the initial density contrast within a 
sphere of radius Xj, centered around a peak in the density field will be 
proportional to the two-point correlation function and will be v 2 £(xi) 
with a probability proportional to exp(— v 2 jl). It follows that the con- 
tribution from a typical region will scale as £ oc while that from 
higher peaks will scale as £ oc £f . In the quasilinear phase, most domi- 
nant contribution arises from high peaks and we find the scaling to be 
Cql oc £f. The non-linear, virialized, regime is dominated by contribu- 
tion from several typical initial regions and has the scaling cjvl oc £^ 2 . 
This was essentially the result obtained in section 2 except that we 
took v = 1. To take into account the statistical fluctuations of the 
initial Gaussian field we can average over different v with a Gaussian 
probability distribution. 

Such an analysis leads to the following result. The relationship between 
£(a, x) and £l(o, V) becomes 



aa,x)=A[Z L (a,l)f h/2 ;A 



20F 



3h/a 



(29) 



where 



6h 



a = — r (30) 

2 + /i(n + 3) V ; 

and A ~ 0.5 is the ratio between the final virialized radius and the 
radius at turn-around. In our model, h = 2 in the quasi-linear regime, 
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and h = 1 in the non-linear regime. However, the above result holds 
for any other value of h. Equation (|29| ) shows that the scaling relations 
(20) acquire coefficients which depend on the spectral index n when we 



average over peaks of different heights. This effect is seen in simulations 
and equation (^) correctly accounts for the numerical results (Munshi 
and Padmanabhan, 1996). 

(ii) In attempting to generalize our results to higher order correlation 
functions, it is important to keep the following aspect in mind. The TVth 
order correlation function will involve N — 1 different length scales. To 
make progress, one needs to assume that, although there are different 
length scales present in reduced n-point correlation function, all of them 
have to be roughly of the same order to give a significant contribution. If 
the correlation functions are described by a single scale, then a natural 
generalisation of equation (g), will be 

Sv - <x?< 7V - 1) )/- 3(iV - 1) (31) 

Given such an ansatz for the TV" point correlation function, one can 
compute the Sn coefficients defined by the relation Sn = £n/%2 1 m 
a straightforward manner. We find that 

S N = (4*r~W \ 2 { (32) 

r(^±i)] 

where a is defined in equation ([3(]). Given the function h(£), this equa- 
tion allows one to compute (approximately) the value of Sn parameters 
in the quasi-linear and non-linear regimes. In our model h = 2 in the 
quasi-linear regime and h = 1 in the non-linear regime. The numerical 
values of Sn computed for different power spectra agrees reasonably 
well with simulation results. (For more details, see Munshi and Pad- 
manabhan, 1996.) 
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